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Abstract. Using lubrication theory, drying processes of sessile colloidal droplets on a solid substrate are 
studied. A simple model is proposed to describe temporal dynamics both the shape of the drop and the 
volume fraction of the colloidal particles inside the drop. The concentration dependence of the viscosity 
is taken into account. It is shown that the final shapes of the drops depend on both the initial volume 
fraction of the colloidal particles and the capillary number. The results of our simulations are in a reasonable 
agreement with the published experimental data. The computations for the drops of aqueous solution of 
human serum albumin (HSA) are presented. 



1 Introduction 

The patterns observed after the complete desiccation of 
sessile drop of a biological fluid such as serum, cerebrospinal 
fluid, urine etc. attract the attention of the researchers at 
least since the 1950s |1I2I3[ . In the last few decades, the 
interest to the pattern formation in drying drops of the 
biological fluids growths stably in the contexts of numer- 
ous applications (e.g. medical tests [41516] . drug screen- 
ing [7], biostabilization [8] etc). State of the art can be 
found in [9TTU] . 

More extensive investigations are performed in the ad- 
jacent field that is evaporative selfassembly of nano- and 
microparticles from an evaporating sessile drop. This pro- 
cess is known as an important surface-patterning tech- 
nique which has potential application in optical, micro- 
electronic, and sensory devices. 

Several models of deposit formation in drying drops 



conservation. The master equation is a diffusion- ad vection 
equation. 

Processes inside desiccated sessile drop depend on var- 
ious physical and chemical conditions such as the wetting 
property of the substrate and surface roughness, the par- 
ticle size, the size distribution, the particle concentration 
in the droplet, the relative humidity of the evaporation 
environment, ionic strength, and pH. 

If the initial volume fraction of the colloidal parti- 
cles inside a drop is large enough, the contact line does 
not recede during the evaporation process, depinning and 
stick-slip motion is not observed. In particular, it was 
shown that desiccation of the colloidal suspensions consid- 
erably depends on the competition between evaporation 
and gelation 24 . Classification of the possible regimes of 
the desiccation can be done in terms of the desiccation 
time td and the gelation time t g . Our investigation is rc- 
stricted to the case when the desiccation time is very short 



*-j were proposed during last few years |ll | 12 | 13 | 14 | 15 | lti | 17 | 18 | ]^g0j^ d with gelation time td « t In this case a solid 



Mostly, the particles were considered as noninteractmg 
and the formed deposit hasn't any effect on both evap- 
oration and bulk flows. It is rather clear that these as- 
sumptions restrict the possible application of the proposed 
models to the particular kind of the solutes. 

One have to utilize quite different approach, if a phase 
transition (sol to gel or sol to glass) occurs as a result of 
increasing the particle volume fraction due to the solvent 
loss and internal flows. The models interpreting deposit 
as impermeable both for hydrodynamical flows and evap- 
oration were proposed recently [22123] . 

Thus, the work [35] is devoted to study the drying 
processes of polymer solutions on a solid substrate en- 
closed by bank. Drying process is studied in the slow limit 
of the solvent evaporation. The model is based on mass 



gelled or glassed foot builds up near the drop edge, while 
the central part remains fluid. 

For this particular situation, the desiccation process 
can be divided into several stages: 

Pregelation stage. The early time regime where a gel-like 
phase does not appear and the whole volume of the 
drop is liquid. During this stage, the solvent evapo- 
rates, the colloidal particles are transferred by the out- 
ward flow to the drop edge and accumulate near the 
contact line. 

Gelation stage. The concentration of colloidal particles at 
the contact line reaches the gelation concentration C g 
above which a gel or glass phase appears. The gel 
(glass) phase region begins to develop from the edge 
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and the phase front monotonically move from edge to 
center. 

Postgelation stage. The drop is solid-like, the loss of sol- 
vent confined among the colloidal particles evaporates 
very slowly, the cracks appeare all around the drop. 

The first stage was investigated in detail in the series 
of works [25 26 27 , as well as in [22;. In the present article, 
we focus efforts on the second stage only. 

The construction of this paper is as follows. In Sec- 
tion^ we describe our model. In Section [3] the obtained 
results are discussed. Finally, we summarize our results 
and conclude this paper. 



2 Model and assumptions 

2.1 Problem Formulation 

We will consider axisymmetric sessile drop of a colloidal 
solution. We assume that volume fraction of the dispersed 
particles is high enough to ensure strong adhesion of the 
drop, in such a way the contact line of a drop is pinned 
during the drying process. When the environmental condi- 
tions are homogeneous and isotropic, the flow field inside 
the drop has to be axisymmetric, too. 

We will use cylindrical coordinates (r, z, <j>), because 
they are natural for the geometry of sessile drops. The 
origin is chosen in the center of the drop on the substrate. 
The coordinate z is normal to the substrate, and the sub- 
strate is described by z = 0, with z being positive on the 
droplet side of the space. The coordinate r is the polar 
radius. Due to the axial symmetry of the problem and our 
choice of the coordinates, no quantity depends on the az- 
imuthal angle. In this case flow field inside a drop has only 
two components: 

v(r, z, t) — u(r, z, t)n r + w(r, z, t)n z . 

The density of particles is supposed to be equal to the 
density of pure solvent. This density, p, is presumed to be 
constant during desiccation. These assumptions are quite 
reasonable for the biological fluids. 

Moreover, we suppose that diffusion transfer is negli- 
gibly compared with advection. 

Only small drops are of the practical interest, thus, we 
assume that the droplet is sufficiently small so that the 
surface tension is dominant, and the gravitational effects 
can be neglected. The temperature gradient along the in- 
terface of the droplet is assumed to be significantly weak 
so that thermal Marangoni flow is not induced. Moreover, 
we will neglect the concentration dependence of the sur- 
face tension a (i.e. the concentration Marangoni flow is 
absent, too). 

The evaporation process can be considered as quasi- 
steady and we will suppose that system is steady at every 
moment. 

Above assumptions allow us to apply the lubrication 
approximation |28j to the drop. This approach was used 
in [T3 22 , too. By applying the lubrication approximation 



to the continuity and Navier- Stokes equations, fluid flow 
is governed by 

dp c 



dr 

1 d(ru) 
r dr 



dz 2 



= 0, 



dz 
di 



dz 



= 0, 



where p is the pressure and rj is the viscosity. 
The boundary conditions should be written as 



P\z=h 



G_d_ 

r dr 



dh 
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dz 



0, 



w\ 



= o, 



(1) 



where h(r,t) is the droplet height profile. 

In the lubrication approximation, the velocity field can 
be written as a function of drop profile h(r,t) 



(if (»))(£_»,), 
rj or \r or \ or J J \ 2 J 

r or \rj or \r or \ or J J \ 6 2 



(2) 



The height-averaged radial velocity calculated using (0) 
is given by 



_ 1 / _ ff_a_d_ fld_ f dh 
h J 3 rj dr \r dr \ dr 



(3) 



From the conservation of solvent, the height-evolution 
equation is given by 



dh _ 1 d(rh {u}) J 
dt r dr p ' 



(4) 



where J(r, t) is the spatial and time-dependent solvent 
mass flux at the air-liquid interface. 

From the conservation of solute, the height-averaged 
solute concentration (C(r,t)) is given by 



d{h(C)) 
dt 



ld_ 

r dr 



(rh(C)(u)). 



(5) 



If a drop is small and surface tension dominates over 
gravity, the shape of the drop is close to an equilibrium 
shape, i.e. spherical cap. Without considerable loss of pre- 
cision we can suppose that a drop has more simple shape 
as a paraboloid of revolution. Nevertheless, the real shape 
of the drop profile near the contact line is questionable. 
Note that widely used assumption about spherical cap 
shape leads to singularities both in evaporation flux and 
flow rate [T3] . Particularly, the precursor film of a wetting 
drop was discussed in [7J. 

The subject of our investigation is only the second 
stage of drying process when colloidal particles are accu- 
mulated at the edge of the droplet. Hence, we can suppose 
that drop edge has very small thickness. 
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Taking into account pinning and the axial symmetry 
of the problem, the boundary and initial conditions for 
the set of equations dU, ([5]) can be written as 



h \t=Q = 
<«>lr=0 = 

dh 
dr 
h 



r=0 
•\r=R = h l 

(C(r,0)} = f(r), 
d(C) 



h { + ho 

(«>lr=fl = 0. 
0. 

h 



(6) 



dr 



= 0, 



where f(r) is the initial concentration of the colloidal par- 
ticle, ho + /if is the initial height of the drop apex, h{ is 
the initial height of the drop edge (/if -C ho). 

To nondimensionalize the equations, the vertical co- 
ordinate was scaled by the initial droplet height ho, the 
radial coordinate was scaled by the drop radius R, the ve- 
locities were scaled by the characteristic viscous velocity 
u c = ^l, where rjo is the viscosity of pure solvent. Time 

was scaled by vapor flux was scaled by J c = 
where k is the thermal conductivity of the liquid, AT is 
the temperature difference between the substrate and sat- 
uration temperatures, L is the heat of vaporization. The 
height-averaged concentration was scaled by gelation con- 
centration C g . 

Upon scaling, the evaporation number, E = and 



where e 



erioL ' 

appear 



the capillary number, Ca = 
as the dimensionless parameters. 

In the dimensionless form, dynamics of the drop profile 
and the concentration of the colloidal particles are gov- 
erned by 

dh _ ld(rh(u)) 



dr 



JE, 



d(h(C)) 



at 



ld_ 

r dr 



(rh(C)(u)), 



where 



(u) 



1 h 2 d /l d 
3Ca ?y dr \r dr 



dh 
dr 



(8) 
(9) 



To simplify notation, we don't use any special symbols for 
the dimensionless quantities. 

In our consideration, we take into account concentra- 
tion dependence of solvent viscosity rj(r, t). Thus, viscosity 
is not a constant but varies in space and time. 

The master equations ([7]) and fl5|) may be reorganized 
and rewritten using the matrix notation 



dF dF 

~dt +{u) ^=^ 



h 



(10) 



(11) 



±^-(r(u))-EJ' 
r dr , , 



(12) 



In the framework of the lubrication approximation, the 
evaporative mass flux must be equal to zero at the con- 
tact line. To our best knowledge the experimental data on 
evaporative flux over the free surface of a colloidal droplet 
with moving phase front inside it are not published yet. 
We will exploit the mass flux of evaporating liquid derived 
from the heat transfer analysis |H] 



J(r,t) = 



1 



h(r,t)' 



where k is a dimensionless nonequilibrium parameter, k — I 
00 for nonvolatile substances and k — > for high volatile 
liquids. A conjecture that colloidal particles accumulated 
at the edge of the droplet must hinder the evaporation at 
the contact line was proposed by Fischer [TJ] . The author 
introduced an additional factor in evaporative flux so that 
the flux decreases exponentially near the contact line. We 
suppose that the density of vapor vanishes not at the edge 
of a drop but at the sol-gel phase boundary. In the other 
words, the vapor flux has to be zero if the concentration 
of the colloidal particles reaches a critical value (C = 1), 
which corresponds to a phase transition (sol to gel or sol 
to glass): 

, . l-(C) 2 (r,t) , , 

J(r,t) = \ ■ 13 

We should notice that the results of computations de- 
pend drastically on particular form of the evaporation 
flux. Our choice (|13[) provides linear decreasing with time 
of the drop mass [12130] and the reasonable dynamics of 
drop shape (see Section 



(7) 2.2 Parameters of the model 



The computations were performed with Ca = 0.01, 0.1, 1, 10 
and E = 0.01,0.1, 1, as well as with the Ca and E, which 
correspond to the drops of biological fluid used in the med- 
ical tests. In the medical tests, the typical drop volume 
varies between 10 and 20 fil with diameter 5-7 mm [SJ 
p. 100]. It means that typical drop height is ho ~ 1 mm. In 
general, the drop height should be lower, since we consider 
only the second stage of desiccation, when the volume frac- 
tion of the colloidal particles reaches the critical value of 
<Pg near the drop edge. Nevertheless, the experiments [31 
as well as the simulations [22 27 show that the first stage 
lasts only about 10 % of desiccation time. Taking into ac- 
count linear decreasing of the drop height, we can ignore 
the height change. We take the middle value R = 3 mm as 
the typical radius, so £ w 0.33. Assuming p = 10 3 kg/m 3 , 
rjo = 10~ 3 Pa-sec, and a = 73 • 10~ 3 N/m, the charac- 
teristic viscous velocity will be u c = 10~ 3 m/sec. In this 
case the capillary number will be Ca w 4 • 10~ 4 . It should 
be noted that e 2 « 0.1, hence the lubrication approxi- 
mation is still valid. The rest of the quantities of use are 
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Fig. 1. Nondimensional viscosity vs. normalized volume 
fraction: experiment [32] and fit by Mooney (HU). 

k = 0.6 W/(m K), AT = 80 K, L = 2.25 ■ 10 6 J/kg, hence 
E w 3 ■ 10~ 2 . 

The experimental data of solution viscosities of HSA 
for concentrations ranges from 8.2 up to 369 kg/m 3 at 
pH= 7.0 and temperature 20° C were taken from [32]. The 
volume fraction of colloidal particles can be calculated as 



</> = 



VC 



where V — ^Trpb 3 is the hydrodynamic volume of one 
dissolved protein (a and b are the effective semi-axes of 
hydrated HSA and p is the axial ratio), Mh = 91.675 kDa 
is the molecular mass of the dissolved protein [32] • 

Taking into account the data published in [33, a = 
8.2 nm, b = 2.1 nm, p = 3.95, we can get the relation 
<P = 10~ 3 C. Then the volume fraction was divided by 
<£ g = 0.5. Nondimensional data were fitted by Mooney's 
formula using least square method 



T) = ?/o exp 



1-K& 



(14) 



where 770 is the viscosity of pure solvent (see fig. H]). 

For the normalized by <P g volume fraction and nondi- 
mensional viscosity, the fitting parameters are S = 2.42 ± 
0.07 and K = 0.84 ± 0.02. 

One can calculate the density of the hydrated protein 
p « 993 kg/m 3 , hence our assumption ([Section l2.ip about 
equality of the solute and solvent densities holds for aque- 
ous solution. 

The initial concentration can be written as 



f(r) = 2 - (Co) + 2 



(Co) - 1 



1 + cxp(w(r — 1)) ' 



where w is an adjustable constant, it is related to the 
length over which the concentration of colloidal particles 
increases rapidly, w — 10 in all computation presented in 
this article. This choice gives us a good agreement with 



0,4 



B 0,3 



Fig. 2. Temporal dynamics of mass, with Ca = 0.1, E 
0.1, <Z> = 0.2. 



the calculated distribution at the end of the first stage 
of desiccation [25126127] . Notice, that /(0) « (C ) and 
/(!) = !• 

Computations in [27] demonstrate, that at the begin 
of the second stage, the concentration in the drop center is 
approximately equal to the initial concentration. Notice, 
that in [22] the concentration in the drop center growths 
at the rate one and a half times. This fact is the natural 
consequence of utilized uniform vapor flux over the free 
surface. 

k = 0.01 and /if = 0.01 in all presented simulations. 
Our computations say that the results are almost insensi- 
tive to /if, if the relation hi <C /io is valid. 



3 Results 

The results of our simulations are presented in figs. [21 El El 
[3 El Z max occurring in the captions relates to the time 
of the whole drop bulk becoming solid-like and the drop 
height almost not decreasing. 

Mass of the drop decreases linearly until almost all 
solvent lost from the drop, then it stays nearly constant 
(Fig. [2]). This behavior agrees with the published exper- 
imental results (see e.g. [12130] ) and demonstrates that 
the proposed evaporation flux (j 13|) is rather adequate. In 
our simulations, difference between mass of the completely 
dried drop and the initial total mass of the colloidal par- 
ticles does not exceed 1.5 %. This fact confirms that com- 
putations are quite correct. 

Drop height decreases approximately in linearly man- 
ner until almost all solvent evaporates from the drop, then 
it remains unchanged (Fig. [3]). This behavior is typical just 
for the case when <§; i g [24 , which is the subject of our 
investigations. 

We have investigated the effect of initial volume frac- 
tion of the colloidal particles inside the drop on time evo- 
lution of the droplet height profile as well as distribution of 
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the colloidal particles (Fig.@}. If the initial volume fraction 
is low there will be visible roller around the edge of dried 
sample. If the initial volume fraction has an intermediate 
value the final shape of the sample looks like a pancake, 
with a depressed central zone. Indeed, if the initial volume 
fraction is close to the value <P g the drop shape varies very 
low. Notice that profile evolution dynamics during early 
time is very similar to the dynamics of the constant base 
model proposed in [11) . 

Fig. [3] shows our results obtained for fixed initial vol- 
ume fraction and evaporation number and with different 
values of capillary number. If capillary number is small 
(i.e. surface-tension forces dominate viscous forces) early 
time evolution of the droplet height profile is similar to 
constant base model [TT]. If the capillary number is large, 
then viscous forces will dominate surface tension forces. 
In this case profile dynamics corresponds to constant an- 
gle model |11) . It looks like if the capillary number is in- 
termediate, then the evolution of profile is closer to the 
predictions of model proposed by Popov Q3]. Thus, the 
known models arc suitable and correct for different par- 
ticular cases and can be reproduced within our model. 

Evolution of the colloidal particle distribution is de- 
picted in Fig. [SI This result is in qualitative agreement 
with 22]. 

Moreover, some computations have been performed for 
the parameters which are close to the values for samples 
used in the medical tests. Time evolution of the droplet 
height profile is presented in Fig. [71 Notice that final shape 
of the drop is a flat pancake, with a depressed central zone 
(FigfTJ) and looks like the real sample (Fig. [5]). 



4 Conclusion 

The main open question of our work is related to the real 
form of evaporation flux and adequacy of our choice (|13p . 
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Fig. 7. (Color online) Time evolution of the droplet height 
profile for t/t max = 0, 0.25, 0.5, 0.75, 1, with Ca = 4 • 10~ 4 , 
E = 0.03, <P = 0.2. 

In our opinion, this question can be resolved only experi- 
mentally. 

The biological fluids contain both colloidal particles 
(proteins) and the salts. Salts have rather high diffusivity, 
hence they cannot be described using ([TUf . It is of interest 
and of value to include in the model additional spices with 
high diffusivity. The natural way to improve the proposed 
model is to supplement the set of equations (jlOD with addi- 
tional advection-diffusion equation describing spatial and 
temporal dynamics of salt [25126127] . 

The authors would like to thank T.A. Yakhno for the photo 
(fig-©- This work was supported by the Russian Foundation 
for Basic Research, project no. 09-08-97010-r_povolzhje_a. 
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a) b) c) 

Fig. 5. (Color online) Time evolution of the droplet height profile for i/t ma x = 0,0.25,0.5,0.75,1, with <P = 0.2, 
E = 0.1. a) Ca = 0.01, b) Ca = 1, c) Ca = 10. 




Fig. 8. Sample of a drop after 16 min it was placed 
onto glass slide (10 w% of BSA in normal saline solution). 
(Courtesy T.A. Yakhno). 
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